Explorar separabilidad

if (!require("sf")) install.packages("sf")
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
if (!require("tidyverse")) install.packages("tidyverse")
## Loading required package: tidyverse
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.2     ✓ dplyr   1.0.6
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
if (!require("GGally")) install.packages("GGally")
## Loading required package: GGally
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
if (!require("np")) install.packages("np")
## Loading required package: np
## Nonparametric Kernel Methods for Mixed Datatypes (version 0.60-10)
## [vignette("np_faq",package="np") provides answers to frequently asked questions]
## [vignette("np",package="np") an overview]
## [vignette("entropy_np",package="np") an overview of entropy-based methods]
source("funiciones_sospechos.R")
e_zona <- st_read("Datos/poligonos_14_16.shp") %>% st_drop_geometry()
## Reading layer `poligonos_14_16' from data source `/home/alequech/Documents/usfsmex/capacitaciones/funcion_densidad_samof/Datos/poligonos_14_16.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1416 features and 59 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 2296366 ymin: 841959 xmax: 2340058 ymax: 872058.3
## Projected CRS: Lambert_Conformal_Conic
names(e_zona)
##  [1] "fid"        "unico"      "OBJECTID"   "OBJECTID_1" "mdmx_t1"   
##  [6] "mdmx_t2"    "ipcc_dyna"  "prob_cmb"   "DNDVI"      "MODO"      
## [11] "SAMOF16"    "Region"     "Reg_Cuadr"  "ID_UTE"     "Val_T1"    
## [16] "Val_T2"     "Metodo"     "Interprete" "OBS"        "Shape_Leng"
## [21] "Shape_Area" "NDVI14"     "NDVI16"     "SUP_MNA"    "Val_Cambio"
## [26] "oid_1"      "change"     "mdmx_t1_2"  "mdmx_t2_2"  "mdmx_t1t2" 
## [31] "mdmx_cls"   "ipcc_dyna_" "mdmx_t1_l"  "prob_cmb_2" "maf_4"     
## [36] "maf_5"      "maf_6"      "maf_1"      "maf_2"      "maf_3"     
## [41] "tcd_t2"     "tcd_t1"     "met_t1_1"   "met_t1_3"   "met_t1_2"  
## [46] "met_t1_5"   "met_t1_4"   "met_t1_7"   "met_t1_6"   "met_t2_6"  
## [51] "met_t2_7"   "met_t2_4"   "met_t2_5"   "met_t2_2"   "met_t2_3"  
## [56] "met_t2_1"   "area"       "id"         "unico_raw"
# tcd tree cover density  * 2015 en los casos posteriores  
# mdmx_t1 etiqueta de cobertura t1 
# mdmx_t2  t2
# ipcc_dyna_ etiqueta 1 a etiqueta 2
# prob_cmb_c prob. de cambio de t1  a t2 basado en la vecindad 
# maf_5 coef. PCA  
# met_t estadicas zonales de bandas 
#t1 = 2014 
#t2 = 2016
# Val_Cambio = C CAMBIO, p 
my_colors <- c("P" = "#00b159","C" = "#d11141")
e_zona %>% 
  select(starts_with("met_") ,Val_Cambio) %>%  
    ggpairs(., aes(colour = Val_Cambio),
        upper = list(continuous = wrap("points", alpha = 0.3)),
        diag = list(discrete="barDiag", continuous = wrap("densityDiag", alpha=0.5 )),
        lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1)),
        legend = c(1,1))+
        scale_color_manual(values = my_colors)+
        scale_fill_manual(values = my_colors)

e_zona %>% 
    filter(SAMOF16 == 12)%>% 
  select(starts_with("met_") ,Val_Cambio) %>% 
    ggpairs(., aes(colour = Val_Cambio),
        upper = list(continuous = wrap("points", alpha = 0.3)),
        diag = list(discrete="barDiag", continuous = wrap("densityDiag", alpha=0.5 )),
        lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1)),
        legend = c(1,1))+
        scale_color_manual(values = my_colors)+
        scale_fill_manual(values = my_colors)

Función de densidad - 1 paso

Todas las categorias

variables <- names(e_zona)[startsWith(names(e_zona),"met_")]
formula <- formula(paste("~", paste(variables, collapse = "+")))
bw <- np::npudensbw(formula, data = e_zona, bwmethod="normal-reference")
f <- npudens(bw)
df_sospechosos  <- cbind(e_zona, f$dens)
df_sospechosos$sospechoso <- with(e_zona, ifelse(f$dens <= quantile(f$dens, 0.10), 1, 0))
ct_t <- table(df_sospechosos$Val_Cambio, df_sospechosos$sospechoso)
ct_t
##    
##        0    1
##   C  108    6
##   P 1165  137
prop.table(ct_t) # cell percentages
##    
##               0           1
##   C 0.076271186 0.004237288
##   P 0.822740113 0.096751412
prop.table(ct_t, 1) # row percentages
##    
##              0          1
##   C 0.94736842 0.05263158
##   P 0.89477727 0.10522273
prop.table(ct_t, 2) # column percentages
##    
##              0          1
##   C 0.08483896 0.04195804
##   P 0.91516104 0.95804196
ct <- table(df_sospechosos$Val_Cambio,df_sospechosos$sospechoso)
ct
##    
##        0    1
##   C  108    6
##   P 1165  137
prop.table(ct) # cell percentages
##    
##               0           1
##   C 0.076271186 0.004237288
##   P 0.822740113 0.096751412
prop.table(ct, 1) # row percentages
##    
##              0          1
##   C 0.94736842 0.05263158
##   P 0.89477727 0.10522273
prop.table(ct, 2) # column percentages
##    
##              0          1
##   C 0.08483896 0.04195804
##   P 0.91516104 0.95804196

para una categoria

sb  <- e_zona %>% 
    filter(SAMOF16 == 12)
bw_sb <- np::npudensbw(formula, data = sb, bwmethod="normal-reference")
f_sb <- npudens(bw_sb)
df_sospechosos_sb  <- cbind(sb, f_sb$dens)
df_sospechosos_sb$sospechoso <- with(sb, ifelse(f_sb$dens <= quantile(f_sb$dens, 0.10), 1, 0))

ct_sd <-  table(df_sospechosos_sb$Val_Cambio,df_sospechosos_sb$sospechoso)
ct_sd
##    
##       0   1
##   C  18   7
##   P 315  30
prop.table(ct_sd, 1)
##    
##              0          1
##   C 0.72000000 0.28000000
##   P 0.91304348 0.08695652
prop.table(ct_sd, 2)
##    
##              0          1
##   C 0.05405405 0.18918919
##   P 0.94594595 0.81081081
# densidad_1paso <-
#   function(df_cat,
#            poutlayer = 0.10,
#            var =  c("navevar1", "navevar2")) {
#     formula <- formula(paste("~", paste(var, collapse = "+")))
#     bw <- np::npudensbw(formula, data = df_cat[,var], bwmethod="normal-reference")
#     #f <- npudens(tdat = df_cat[, var])
#     f <- np::npudens(bw)
#     df_cat$fdens <- f$dens
#     df_cat$sospechoso <- with(df_cat, ifelse(df_cat$fdens <= quantile(df_cat$fdens, poutlayer), 1, 0))
#     return(df_cat)
#   }
var_t1 <-  names(sb)[startsWith(names(sb),"met_t1")]
sb_t1 <- densidad_1paso(sb ,poutlayer = 0.10, var = var_t1) 
ct_sd_1 <-  table(sb_t1$Val_Cambio,sb_t1$sospechoso)
ct_sd_1
##    
##       0   1
##   C  20   5
##   P 313  32
prop.table(ct_sd_1, 1)
##    
##              0          1
##   C 0.80000000 0.20000000
##   P 0.90724638 0.09275362
prop.table(ct_sd_1, 2)
##    
##              0          1
##   C 0.06006006 0.13513514
##   P 0.93993994 0.86486486

Función de densidad multiples iteraciones

sb_mp  <- densidad_m_pasos(df_cat = sb ,poutlayer = 0.10, var = variables, numit = 10) 
ct_sd_mp <-  table(sb_mp$Val_Cambio,sb_mp$sospechoso)
ct_sd_mp
##    
##       0   1
##   C  18   7
##   P 312  33
prop.table(ct_sd_mp, 1)
##    
##              0          1
##   C 0.72000000 0.28000000
##   P 0.90434783 0.09565217
prop.table(ct_sd_mp, 2)
##    
##              0          1
##   C 0.05454545 0.17500000
##   P 0.94545455 0.82500000